#include "filter.h"
#include <iostream>
using namespace std;


void Filter::setKernel(double *kernel,int k_rows,int k_cols){
    this->kernel=kernel;
    this->k_rows=k_rows;
    this->k_cols=k_cols;
//    double divisor=0;
//    for(int i=0;i<k_cols;i++)
//        for(int j=0;j<k_rows;j++)
//            divisor+=kernel[i+j*k_cols];
//    for(int i=0;i<k_cols;i++)
//        for(int j=0;j<k_rows;j++)
//            kernel[i+j*k_cols]/=divisor;
}

void Filter::setKernel(double *kernel,int k_rows,int k_cols, double divisor){
    this->kernel=kernel;
    this->k_rows=k_rows;
    this->k_cols=k_cols;
    for(int i=0;i<k_cols;i++)
        for(int j=0;j<k_rows;j++)
            kernel[i+j*k_cols]/=divisor;
}

QImage Filter::normalize(int *R,int *G,int *B){

    return normalize(R,G,B,0);
    }

QImage Filter::normalize(int* R,int* G,int* B,int mode) //-1 cut negatives, 0 let them be, 1 abs of neg.
{
    return normalize(R,G,B,mode,0.0);
}
QImage Filter::normalize(int* R,int mode,double frac) //-1 cut negatives, 0 let them be, 1 abs of neg.
{                                                     //-2 no normalization
    int rows=in_image->height();
    int cols=in_image->width();
    if(mode==-2)
    {
        QImage nonorm(in_image->width(),in_image->height(),in_image->format());
        for(int i=0;i<cols;i++)
            for(int j=0;j<rows;j++)
            {
                int r=R[i+cols*j];
                if(r<0) r=0;
                if(r>255) r=255;
                nonorm.setPixel(i,j,qRgb(r,r,r));
            }
         return nonorm;
    }
    else if(mode==-1)
    for(int i=0;i<cols;i++)
        for(int j=0;j<rows;j++)
        {
            if(R[i+j*cols]<0)
                R[i+j*cols]=0;
        }
    else if(mode==1)
        for(int i=0;i<cols;i++)
            for(int j=0;j<rows;j++)
            {
                if(R[i+j*cols]<0)
                    R[i+j*cols]=abs(R[i+j*cols]);
            }
    int max=R[0],min=R[0];
    for(int i=0;i<cols;i++)
        for(int j=0;j<rows;j++){
        if(R[i+j*cols]>max)
            max=R[i+j*cols];
        if(R[i+j*cols]<min)
            min=R[i+j*cols];
    }
    int bias=0;
    if(min<0) bias=abs(min);
    int *histR=new int[max+1+bias];
    for(int i=0;i<max+1+bias;i++)
    {
        histR[i]=0;
    }
    for(int i=0; i<cols;i++)
        for(int j=0;j<rows;j++)
        {
            histR[R[i+j*cols]+bias]++;
        }
    int maxhistR=0;
    for(int i=0;i<max+1+bias;i++)
    {
        if(histR[i]>maxhistR)
            maxhistR=histR[i];
    }
    int fraccR=frac*maxhistR;
    for(int i=0;i<max+1+bias;i++)
        if(histR[i]>fraccR)
        {
            min=i-bias;
            break;
        }
    for(int i=max+bias;i>=0;i--)
        if(histR[i]>fraccR)
        {
            max=i-bias;
            break;
        }
    QImage qout=QImage(in_image->width(),in_image->height(),in_image->format());
    double m=255.0/(max-min);
    double c=-m*min;
    cout<<"m is"<<m<<"max is"<<max<<"min is"<<min<<endl;
    for(int i=0;i<cols;i++)
        for(int j=0;j<rows;j++){
        int r=m*R[i+j*cols]+c;
        if(r<0) r=0;
        if(r>255) r=255;
        qout.setPixel(i,j,qRgb(r,r,r));
    }
    delete [] R;
    delete []histR;
    return qout;


}
QImage Filter::normalize(int* R,int* G,int* B,int mode,double frac) //-1 cut negatives, 0 let them be, 1 abs of neg.
{                                                                   //-2 no normalization
    int rows=in_image->height();
    int cols=in_image->width();
    if(mode==-2)
    {
        QImage nonorm(in_image->width(),in_image->height(),in_image->format());
        for(int i=0;i<cols;i++)
            for(int j=0;j<rows;j++)
            {
                int r=R[i+cols*j];
                if(r<0) r=0;
                if(r>255) r=255;
                int g=G[i+cols*j];
                if(g<0) g=0;
                if(g>255) g=255;
                int b=B[i+cols*j];
                if(b<0) b=0;
                if(b>255) b=255;
                nonorm.setPixel(i,j,qRgb(r,g,b));
            }
         return nonorm;
    }
    else if(mode==-1)
    for(int i=0;i<cols;i++)
        for(int j=0;j<rows;j++)
        {
            if(R[i+j*cols]<0)
                R[i+j*cols]=0;
            if(G[i+j*cols]<0)
                G[i+j*cols]=0;
            if(B[i+j*cols]<0)
                B[i+j*cols]=0;
        }
    else if(mode==1)
        for(int i=0;i<cols;i++)
            for(int j=0;j<rows;j++)
            {
                if(R[i+j*cols]<0)
                    R[i+j*cols]=abs(R[i+j*cols]);
                if(G[i+j*cols]<0)
                    G[i+j*cols]=abs(G[i+j*cols]);
                if(B[i+j*cols]<0)
                    B[i+j*cols]=abs(B[i+j*cols]);
            }
    int max=R[0],min=R[0];
    for(int i=0;i<cols;i++)
        for(int j=0;j<rows;j++){
        if(R[i+j*cols]>max)
            max=R[i+j*cols];
        if(G[i+j*cols]>max)
            max=G[i+j*cols];
        if(B[i+j*cols]>max)
            max=B[i+j*cols];
        if(R[i+j*cols]<min)
            min=R[i+j*cols];
        if(G[i+j*cols]<min)
            min=G[i+j*cols];
        if(B[i+j*cols]<min)
            min=B[i+j*cols];
    }
    int bias=0;
    if(min<0) bias=abs(min);
    int *histR=new int[max+1+bias];
    int *histG=new int[max+1+bias];
    int *histB=new int[max+1+bias];
    for(int i=0;i<max+1+bias;i++)
    {
        histR[i]=0;
        histG[i]=0;
        histB[i]=0;
    }
    for(int i=0; i<cols;i++)
        for(int j=0;j<rows;j++)
        {
            histR[R[i+j*cols]+bias]++;
            histG[G[i+j*cols]+bias]++;
            histB[B[i+j*cols]+bias]++;
        }
    int maxhistR=0;
    //int maxhistG=0;
    //int maxhistB=0;
    for(int i=0;i<max+1+bias;i++)
    {
        if(histR[i]>maxhistR)
            maxhistR=histR[i];
        if(histG[i]>maxhistR)
            maxhistR=histG[i];
        if(histB[i]>maxhistR)
            maxhistR=histB[i];
    }
    int fraccR=frac*maxhistR;
    //int fraccG=fracc*maxhistG;
    //int fraccB=fracc*maxhistB;
    for(int i=0;i<max+1+bias;i++)
        if(histR[i]>fraccR || histG[i]>fraccR || histB[i]>fraccR)
        {
            min=i-bias;
            break;
        }
    for(int i=max+bias;i>=0;i--)
        if(histR[i]>fraccR || histG[i]>fraccR || histB[i]>fraccR)
        {
            max=i-bias;
            break;
        }
    QImage qout=QImage(in_image->width(),in_image->height(),in_image->format());
    double m=255.0/(max-min);
    double c=-m*min;
    cout<<"m is"<<m<<"max is"<<max<<"min is"<<min<<endl;
    for(int i=0;i<cols;i++)
        for(int j=0;j<rows;j++){
        int r=m*R[i+j*cols]+c;
        if(r<0) r=0;
        if(r>255) r=255;
        int g=m*G[i+j*cols]+c;
        if(g<0) g=0;
        if(g>255) g=255;
        int b=m*B[i+j*cols]+c;
        if(b<0) b=0;
        if(b>255) b=255;
        qout.setPixel(i,j,qRgb(r,g,b));
    }
    delete [] R;
    delete []histR;
    delete []histG;
    delete []histB;
    delete [] G;
    delete [] B;
    return qout;


}


struct image Filter::convolve_img(QImage *in_image){
    int i,j,a,b,l,k;
    a=this->k_rows/2;
    b=this->k_cols/2;
    this->in_image=in_image;
    int cols=in_image->width();
    int rows=in_image->height();
    int *out_arrayR=new int[cols*rows];
    int *out_arrayG=new int[cols*rows];
    int *out_arrayB=new int[cols*rows];
    for(j=0;j<rows;j++)
    for(i=0;i<cols;i++)
    {
          out_arrayR[i+j*cols]=0;
          out_arrayG[i+j*cols]=0;
          out_arrayB[i+j*cols]=0;
          for(l=0;l<k_rows;l++)
          for(k=0;k<k_cols;k++)
          {
                  int jj=j+l-a;
                  int ii=i+k-b;
                  if(jj>=0 && ii>=0 && jj<in_image->height() && ii<in_image->width()){
                      int kernel_pos=this->k_rows-k-1+this->k_cols*(this->k_rows-l-1);
                      double ker_val=*(this->kernel+kernel_pos);
                      out_arrayR[i+j*cols]=out_arrayR[i+j*cols]+ker_val*qRed(in_image->pixel(ii,jj));
                      out_arrayG[i+j*cols]=out_arrayG[i+j*cols]+ker_val*qGreen(in_image->pixel(ii,jj));
                      out_arrayB[i+j*cols]=out_arrayB[i+j*cols]+ker_val*qBlue(in_image->pixel(ii,jj));
                  }
          }
    }
  struct image img;
  img.R=out_arrayR;
  img.G=out_arrayG;
  img.B=out_arrayB;
  return img;
}

QImage Filter::convolve(QImage *in_image, int mode)
{
   return convolve(in_image,mode,0);
}

QImage Filter::convolve(QImage* in_image,int mode, double frac){
    this->in_image=in_image;
    int cols=in_image->width();
    int rows=in_image->height();
    int *inputR=new int[cols*rows];
    int *inputG;
    int *inputB;
    int gray=in_image->isGrayscale();
    if(!(gray))
    {
        inputG=new int[cols*rows];
        inputB=new int[cols*rows];
    }
    for(int i=0;i<cols;i++)
        for(int j=0;j<rows;j++)
        {
            inputR[i+j*cols]=qRed(in_image->pixel(i,j));
            if(!(gray))
            {
            inputG[i+j*cols]=qGreen(in_image->pixel(i,j));
            inputB[i+j*cols]=qBlue(in_image->pixel(i,j));
            }
        }
    int* outputR=convolve(inputR);
    delete []inputR;
    int* outputG;
    int* outputB;
    if(!(gray))
    {
        outputG=convolve(inputG);
        delete []inputG;
        outputB=convolve(inputB);
        delete []inputB;
    }
    if(gray)
    {
        return normalize(outputR,mode,frac);
    }
  return normalize(outputR,outputG,outputB,mode,frac);
}

QImage Filter::convolve(QImage *in_image){

    this->in_image=in_image;
    int cols=in_image->width();
    int rows=in_image->height();
    int *inputR=new int[cols*rows];
    int *inputG;
    int *inputB;
    int gray=in_image->isGrayscale();
    if(!gray)
    {
        inputG=new int[cols*rows];
        inputB=new int[cols*rows];
    }
    for(int i=0;i<cols;i++)
        for(int j=0;j<rows;j++)
        {
            inputR[i+j*cols]=qRed(in_image->pixel(i,j));
            if(!(gray))
            {
            inputG[i+j*cols]=qGreen(in_image->pixel(i,j));
            inputB[i+j*cols]=qBlue(in_image->pixel(i,j));
            }
        }
    int* outputR=convolve(inputR);
    delete []inputR;
    int* outputG;
    int* outputB;
    if(!(gray))
    {
        outputG=convolve(inputG);
        delete []inputG;
        outputB=convolve(inputB);
        delete []inputB;
    }
  QImage newimage(in_image->width(),in_image->height(),in_image->format());
  for(int i=0;i<cols;i++)
          for(int j=0;j<rows;j++)
          {
              int r=outputR[i+cols*j];
              if(r<0) r=0;
              if(r>255) r=255;
              if(!(gray))
              {
              int g=outputG[i+cols*j];
              if(g<0) g=0;
              if(g>255) g=255;
              int b=outputB[i+cols*j];
              if(b<0) b=0;
              if(b>255) b=255;
              newimage.setPixel(i,j,qRgb(r,g,b));
              }else
              newimage.setPixel(i,j,qRgb(r,r,r));
          }
  delete []outputR;
  if(!(gray))
  {
  delete []outputG;
  delete []outputB;
  }
  return newimage;
}

int * Filter::convolve(int * input)
{
    int a=this->k_rows/2;
    int b=this->k_cols/2;
    int cols=in_image->width();
    int rows=in_image->height();
    int *output=new int[cols*rows];
    for(int j=0;j<rows;j++)
    for(int i=0;i<cols;i++)
    {
          double sum=0;
          output[i+j*cols]=0;
          for(int l=0;l<k_rows;l++)
          for(int k=0;k<k_cols;k++)
          {
                  int jj=j+l-a;
                  int ii=i+k-b;
                  if(jj<0) jj=0;
                  if(jj>=in_image->height()) jj=in_image->height()-1;
                  if(ii<0) ii=0;
                  if(ii>=in_image->width()) ii=in_image->width()-1;
                  int kernel_pos=this->k_rows-k-1+this->k_cols*(this->k_rows-l-1);
                  double ker_val=*(this->kernel+kernel_pos);
                  sum+=ker_val*input[ii+jj*cols];

          }
          output[i+j*cols]=(int)sum;
    }
    return output;

}

QImage Filter::normalize(QImage *qin){
QImage qout=QImage(qin->width(),qin->height(),qin->format());
int rows=qin->height();
int cols=qin->width();
int i,j,max=0,min=255,R,G,B;
for(i=0;i<cols;i++)
    for(j=0;j<rows;j++){
    R=qRed(qin->pixel(i,j));
    G=qGreen(qin->pixel(i,j));
    B=qBlue(qin->pixel(i,j));
    if(R>max)
        max=R;
    if(G>max)
        max=G;
    if(B>max)
        max=B;
    if(R<min)
        min=R;
    if(G<min)
        min=G;
    if(B<min)
        min=B;
}
double m=255.0/(max-min);
double c=-m*min;
cout<<"m is"<<m<<"max is"<<max<<"min is"<<min<<endl;
for(i=0;i<cols;i++)
    for(j=0;j<rows;j++){
    R=qRed(qin->pixel(i,j));
    G=qGreen(qin->pixel(i,j));
    B=qBlue(qin->pixel(i,j));
    qout.setPixel(i,j,qRgb(m*R+c,m*G+c,m*B+c));
}
return qout;
}

void Filter::fft(double *real, double *imag, int size, int dir)
{
    int npoints;
   // Number of points needed
   int times=0;
   for (npoints=1;npoints<size;times++)
      npoints *= 2;

   // Reversal
   int i2 = npoints >> 1;
   int j = 0;
   int k;
   for (int i=0;i<npoints-1;i++) {
      if (i < j) {
         double x = real[i];
         double y = imag[i];
         real[i] = real[j];
         imag[i] = imag[j];
         real[j] = x;
         imag[j] = y;
      }
      k = i2;
      while (k <= j) {
         j -= k;
         k >>= 1;
      }
      j += k;
   }

   //Compute the FFT
   double c1 = -1.0;
   double c2 = 0.0;
   int l2 = 1;
   for (int l=0;l<times;l++) {
      int l1 = l2;
      l2 <<= 1;
      double u1 = 1.0;
      double u2 = 0.0;
      for (int j=0;j<l1;j++) {
         for (int i=j;i<npoints;i+=l2) {
            int i1 = i + l1;
            double t1 = u1 * real[i1] - u2 * imag[i1];
            double t2 = u1 * imag[i1] + u2 * real[i1];
            real[i1] = real[i] - t1;
            imag[i1] = imag[i] - t2;
            real[i] += t1;
            imag[i] += t2;
         }
         double z =  u1 * c1 - u2 * c2;
         u2 = u1 * c2 + u2 * c1;
         u1 = z;
      }
      c2 = sqrt((1.0 - c1) / 2.0);
      if (dir == 1)
         c2 = -c2;
      c1 = sqrt((1.0 + c1) / 2.0);
   }

   //Scaling on the inverse
   if (dir == -1) {
      for (int i=0;i<npoints;i++) {
         real[i] = real[i]/npoints;
         imag[i] = imag[i]/npoints;
      }
   }
}

void Filter::fft2D(double **re, double** im,int width,int height,int dir)
{
   double *real,*imag;

   //Rows
   real = new double[width * sizeof(double)];
   imag = new double[width * sizeof(double)];

   for (int j=0;j<height;j++) {
      for (int i=0;i<width;i++) {
         real[i] = re[i][j];
         imag[i] = im[i][j];
      }
      fft(real,imag,width,dir);
      for (int i=0;i<width;i++) {
         re[i][j] = real[i];
         im[i][j] = imag[i];
      }
   }
   delete [] real;
   delete [] imag;

   //Columns
   real = new double[height * sizeof(double)];
   imag = new double[height * sizeof(double)];

   for (int i=0;i<width;i++) {
      for (int j=0;j<height;j++) {
         real[j] = re[i][j];
         imag[j] = im[i][j];
      }
      fft(real,imag,height,dir);
      for (int j=0;j<height;j++) {
         re[i][j] = real[j];
         im[i][j] = imag[j];
      }
   }
   delete [] real;
   delete [] imag;
}

void Filter::fft2D(int **re, int** im,int width,int height,int dir)
{
   double *real,*imag;

   //Rows
   real = new double[width * sizeof(double)];
   imag = new double[width * sizeof(double)];

   for (int j=0;j<height;j++) {
      for (int i=0;i<width;i++) {
         real[i] = re[i][j];
         imag[i] = im[i][j];
      }
      fft(real,imag,width,dir);
      for (int i=0;i<width;i++) {
         re[i][j] = real[i];
         im[i][j] = imag[i];
      }
   }
   delete [] real;
   delete [] imag;

   //Columns
   real = new double[height * sizeof(double)];
   imag = new double[height * sizeof(double)];

   for (int i=0;i<width;i++) {
      for (int j=0;j<height;j++) {
         real[j] = re[i][j];
         imag[j] = im[i][j];
      }
      fft(real,imag,height,dir);
      for (int j=0;j<height;j++) {
         re[i][j] = real[j];
         im[i][j] = imag[j];
      }
   }
   delete [] real;
   delete [] imag;
}

QImage Filter::convolveFFT(QImage* input)
{
    in_image=input;
    double **kernel=new double*[k_cols];

    for(int i=0;i<k_cols;i++)
    {
        kernel[i]=new double[k_rows];
        for(int j=0;j<k_rows;j++)
            kernel[i][j]=this->kernel[i+k_cols*j];
    }

    int m,n,m2,n2;
    for(m=1;m<in_image->width();m*=2) ;
    for(n=1;n<in_image->height();n*=2) ;
    for(m2=1;m2<k_cols;m2*=2) ;
    for(n2=1;n2<k_rows;n2*=2) ;
    if(m2>m) m=m2;
    if(n2>n) n=n2;
    double ** rimage=new double*[m];
    double ** iimage=new double*[m];
    double ** rkernel=new double*[m];
    double ** ikernel=new double*[m];

    for(int i=0;i<m;i++)
    {
        rimage[i]=new double[n];
        iimage[i]=new double[n];
        rkernel[i]=new double[n];
        ikernel[i]=new double[n];
    }

    //RED
    for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
        {
            if(i<in_image->width() && j<in_image->height())
                rimage[i][j]=(double)qRed(in_image->pixel(i,j));
            else
                rimage[i][j]=0.0;
            iimage[i][j]=0.0;
            if(i<k_cols && j<k_rows)
                rkernel[i][j]=kernel[i][j];
            else
                rkernel[i][j]=0.0;
            ikernel[i][j]=0.0;
        }

    fft2D(rimage,iimage,m,n,1);
    fft2D(rkernel,ikernel,m,n,1);
    for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
        {

            int rval=rimage[i][j]*rkernel[i][j]-iimage[i][j]*ikernel[i][j];
            int ival=rimage[i][j]*ikernel[i][j]+iimage[i][j]*rkernel[i][j];
            rimage[i][j]=rval;
            iimage[i][j]=ival;
        }

    fft2D(rimage,iimage,m,n,-1);
    int *r=new int[in_image->width()*in_image->height()];
    if(in_image->isGrayscale())
    {

        for(int i=0;i<in_image->width();i++)
            for(int j=0;j<in_image->height();j++)
            {
                r[i+in_image->width()*j]=(int)rimage[i][j];
            }
        for(int i=0; i<m; i++)
        {
            delete [] rimage[i];
            delete [] iimage[i];
            delete [] rkernel[i];
            delete [] ikernel[i];
        }
        delete rimage;
        delete iimage;
        delete rkernel;
        delete ikernel;
        return normalize(r,-2,0.0);
    }
    int *g=new int[in_image->width()*in_image->height()];
    int *b=new int[in_image->width()*in_image->height()];
    for(int i=0;i<in_image->width();i++)
        for(int j=0;j<in_image->height();j++)
            r[i+in_image->width()*j]=(int)rimage[i][j];
    //GREEN
    for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
        {
            if(i<in_image->width() && j<in_image->height())
                rimage[i][j]=(double)qGreen(in_image->pixel(i,j));
            else
                rimage[i][j]=0.0;
            iimage[i][j]=0.0;
            if(i<k_cols && j<k_rows)
                rkernel[i][j]=kernel[i][j];
            else
                rkernel[i][j]=0.0;
            ikernel[i][j]=0.0;
        }

    fft2D(rimage,iimage,m,n,1);
    fft2D(rkernel,ikernel,m,n,1);
    for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
        {

            int rval=rimage[i][j]*rkernel[i][j]-iimage[i][j]*ikernel[i][j];
            int ival=rimage[i][j]*ikernel[i][j]+iimage[i][j]*rkernel[i][j];
            rimage[i][j]=rval;
            iimage[i][j]=ival;
        }

    fft2D(rimage,iimage,m,n,-1);
    for(int i=0;i<in_image->width();i++)
        for(int j=0;j<in_image->height();j++)
            g[i+in_image->width()*j]=(int)rimage[i][j];


    //BLUE
    for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
        {
            if(i<in_image->width() && j<in_image->height())
                rimage[i][j]=(double)qBlue(in_image->pixel(i,j));
            else
                rimage[i][j]=0.0;
            iimage[i][j]=0.0;
            if(i<k_cols && j<k_rows)
                rkernel[i][j]=kernel[i][j];
            else
                rkernel[i][j]=0.0;
            ikernel[i][j]=0.0;
        }

    fft2D(rimage,iimage,m,n,1);
    fft2D(rkernel,ikernel,m,n,1);
    for(int i=0;i<m;i++)
        for(int j=0;j<n;j++)
        {

            int rval=rimage[i][j]*rkernel[i][j]-iimage[i][j]*ikernel[i][j];
            int ival=rimage[i][j]*ikernel[i][j]+iimage[i][j]*rkernel[i][j];
            rimage[i][j]=rval;
            iimage[i][j]=ival;
        }

    fft2D(rimage,iimage,m,n,-1);
    for(int i=0;i<in_image->width();i++)
        for(int j=0;j<in_image->height();j++)
            b[i+in_image->width()*j]=(int)rimage[i][j];
    for(int i=0; i<m; i++)
    {
        delete [] rimage[i];
        delete [] iimage[i];
        delete [] rkernel[i];
        delete [] ikernel[i];
    }
    delete rimage;
    delete iimage;
    delete rkernel;
    delete ikernel;
    return normalize(r,g,b,-2,0.0);
}


Filter::~Filter(){


}
